Mott transition of fermionic atoms in a three-dimensional optical trap 
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We study theoretically the Mott metal-insulator transition for a system of fermionic atoms con- 
fined in a three-dimensional optical lattice and a harmonic trap. We describe an inhomogeneous 
system of several thousand sites using an adaptation of dynamical mean field theory solved efficiently 
with the numerical renormalization group method. Above a critical value of the on-site interaction, 
a Mott-insulating phase appears in the system. We investigate signatures of the Mott phase in the 
density profile and in time-of-fiight experiments. 
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Introduction: Ultracold atoms in optical lattices offer 
exciting possibilities to investigate many-particle effects 
and to realize and measure models of condensed mat- 
ter physics like the Hubbard model with unprecedented 
control of the band structure and interaction strength. 

One of the most dramatic effects of strong correlations 
is the Mott transition, where strong interactions drive a 
system insulating. While for bosonic atoms j,] the Mott 
transition in an optical lattice has been realized a few 
years ago by Greiner et al. [3] , the corresponding exper- 
iment for fermionic systems turns out to be much more 
difficult, partially due to problems with cooling and due 
to the need to work with two fermionic species to model 
the spin degree of freedom. However, due to the enor- 
mous experimental progress, a realization of the Mott 
transition in fermionic systems is expected in the near 
future. For example, recently Kohl et al. [3] succeeded 
to capture a single species of fcrmions in an optical lat- 
tice and to observe the Fermi surface. Also two species 
of fermions have been successfully trapped and cooled 
down in a series of experiments studying e.g. the BEC- 
BCS crossover Q. 

For the interpretation of the experiments it is essential 
to investigate the effects of the smooth external confining 
potential holding the atoms in the trap. The resulting 
inhomogeneities can make it more difficult to interpret 
the experiments, but induce also new interesting effects, 
e.g. associated to the sharp surface between metallic 
and insulating regions. A frequently used approximation 
(e.g. in the present context by Ref. is to describe the 
trapped atoms locally by a homogeneous system as in the 
local density approximation. However, in the presence of 
sharp domain walls between two phases, such an approx- 
imation is expected to fail and a more realistic treatment 
of the inhomogeneous system is necessary. 

For one-dimensional (ID) systems, powerful numerical 
d, 01 and analytical Q methods exists to study theoret- 
ically the Mott transition in fermionic systems. For ex- 
ample, in Refs. [5i, 6, 7] quantum Monte Carlo techniques 
were used to investigate the signatures of Mott phases in 
the presence of an external harmonic confinement poten- 



tial for a ID system. Rigol et al. [1,01 argued, that in one 
dimension the inhomogeneities resulting from the trap- 
ping potential essentially destroy the main signatures of 
Mott phases in time-of-flight experiments. 

As exact numerical methods for fermions can only be 
applied to very small systems, one has to resort to ap- 
proximations to calculate the properties of three dimen- 
sional lattices of realistic size. Here the method of choice 
is the so-called dynamical mean field theory (DMFT) 
[1, 0|- Within DMFT, the only approximation is to ne- 
glect non-local contributions to the self-energy. This al- 
lows to map the N site lattice problems to N single- 
impurity Anderson models coupled by a self-consistency 
condition, see discussion below. DMFT is, for example, 
frequently used to describe complex bulk materials, e.g. 
by combining DMFT with band-structure calculations to 
obtain an ab-inito description of strongly correlated ma- 
terials. In a few cases, DMFT has been employed to 
describe inhomogeneous systems H, 11, 12, 13 1 like the 



111 ] or disordered materials 



surface of Mott insulators 

M- 

A main problem of DMFT is the need for an reliable 
and efficient method to solve the effective impurity prob- 
lems. Previous applications of DMFT to inhomogeneous 
systems were using impurity solvers like a two-site ap- 



proximation ll| or slave-boson mean-field theory 14| . 
implying severe further approximations, or started from 
simplified fermionic models such as the Falicov-Kimball 
model Ts^ . We will show that one can also use efficiently 
one of the most accurate irnmirity solvers, the numerical 
renormalization group [illli (NRG), to obtain rehable 



results for traps containing several thousand atoms mod- 
eling a fermionic Hubbard model. 

After introducing the model and our method (DMFT 
for inhomogeneous systems-|-NRG), we will show the re- 
sulting spectral functions and discuss how the transition 
from a metal to a Mott insulating phase can be seen in 
real-space and time-of-flight experiments. We investigate 
the role of temperature, filling and interaction strength. 

Model and Method: We consider the fermionic Hub- 
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FIG. 1; (color online) The number of fermions (n^) per lattice 
site, as a function of the distance r to the origin. The signa- 
ture of the Mott insulating phase is the presence of a plateau 
with (ui) = 1. Note that there are inequivalent sites with dif- 
ferent occupation but the same r. Main panel: Crossover from 
weak to strong interactions for a fixed number oi N = 2869 
particles in the trap at T = 0, [/ = 6J. Upper inset: Depen- 
dence on the number of particles for U = 4.5-D, T — 0. Lower 
inset: T dependence for U = 4.5D, iV = 2869 (the T = and 
T = 0.13-D curves lie on top of each other). 
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FIG. 2: (color online) Local spectral functions in the Mott 
insulating phase {U = 4.5D, iV = 2869) for two different T 
for lattice points with coordinates (x, 1/2, 1/2). Left inset: 
The coherence peak at the Fermi energy, characteristic for a 
strongly correlated metal, vanishes with increasing T. Right 
inset: Close to edge of the atomic cloud, where the potential 
becomes steep, finite size effects are visible. 



bard model on a 3D-cubic lattice, 



u 



+ Vo^rln,,, (1) 



where q^. creates a fermion at site i with spin a, rii^ — 
cl^Cia is the local density, J the nearest-neighbor tunnel- 
ing matrix element and U the effective on-site interac- 
tion and the lattice distance is set to unity. We include 
4224 sites with a distance ri < R — 10 from the center 
of the trap, located in the middle of 8 central sites. The 
strength of the confining harmonic potential Vq — 0.276 J 
is chosen such that all sites with ri > R are unoccupied 
and can be neglected. Using the symmetry of the cubic 
lattice, one has to deal with only 118 inequivalent sites. 

The basic idea of DMFT ^ is to pick a single site 
of the lattice (the 'impurity' i) and model the effect of 
all the other sites by a non-interacting bath of fermions. 
The resulting Anderson impurity model is defined by the 
local interaction U and the hybridization of the impurity 
to the bath. The latter is encoded in the U — Green's 
function, G^ndi(i^)- For this model, one determines the 
local self-energy ^i{uj) (see below). From the Ei(aj), one 
can construct the lattice Green's function 

Giat~^(t^)y = <5ij (cj + Ai - Si(a;) - Var^) - , (2) 

where Jij = J if sites i and j are nearest neighbors and 
otherwise. The bath of each impurity is then determined 
from the requirement that at each site the lattice Green's 
function and the Green's function of the impurity model 



coincide, Giat(w)ii = t/And,i('^) ^ - Sj(w) , thus es- 
tablishing a self-consistency loop. The scheme described 
above scheme can be derived using as the only approxi- 
mation that the self-energy is a local quantity. Both the 
non-local single particle quantum mechanics of fermions 
and all local effects of strong interactions are correctly 
described within DMFT. 

A main difficulty of DMFT is, however, an accurate 
calculation of the self-energy of the Anderson impurity 



model. For this we use the NRG [18|, see Ref. [16| for 
a description of the method. To obtain efficiently Giat 
from an inversion of Eq. ^ it is essential to use the full 
symmetry of the cubic lattice. 

In this paper we restrict ourselves to paramagnetic so- 
lutions which simplifies the rather challenging numerics 
considerably. Also experimentally, it is very difficult to 
reach the low temperatures below which magnetism is 
expected. Furthermore, we do not expect that magnetic 
order will change the density profiles or time-of-flight pic- 
tures considerably. 

Results: Fig. [T] shows how the number of fermions 
per site evolves for increasing interactions U/J, which 
push the fermions away from the center of the trap. 
For the chosen parameters, we obtain for U/D — 0, 1, 2 
{D — 6J is half the bandwidth) a band insulator in 
the center of the trap and a metal further outside. For 
the homogeneous system the critical interaction is given 
by Uc/D = 2.52, but already for U/D = 2 the com- 
pressibility close to half-filling is strongly reduced as can 
be seen in a shoulder in the curve for (n^) w 1. For 
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FIG. 3: (color online) Momentum distribution n^°^ for U / D = 
0, 2, 3, 4.5. Both in the predominantly band-insulating phase 
[U/D = 0) and Mott insulating phase {U/D = 4.5) the 
curves are considerably flatter than for U / D = 2,3 where 
most fermions are in the metallic phase. 



U/D > 3 the incompressibility of the Mott insulating 
state, dn/dfi — 0, manifests itself in a plateau. The 
thickness of this Mott insulating 'onion shell' increases 
for increasing U eliminating the metallic phase in the 
center for U/D = 4.5. The insets of Fig. [T] show how the 
Mott insulating region aXU /D = 4.5 shrinks again when 
the number of fermions in the trap is reduced and how 
thermal excitations destroy the Mott plateau. 

Our method allows to study the spatial dependence of 
the spectral functions. While this quantity is difficult 
to measure for atoms in a trap, it is a highly sensitive 
probe of the metal- insulator transition. Fig [5] shows how 
the local spectral function evolves when moving from the 
center to the edge of the trap at U/D = 4.5. In the 
insulating regime, the spectral function A{uj) is charac- 
terized by the two Hubbard bands with equal weight and 
A{lij = 0) becomes very small (it never vanishes exactly 
as atoms from the metallic regions can tunnel into the 
insulator). As a function of the distance from the center, 
the Hubbard bands shift due to the harmonic potential. 
When the Fermi energy starts to merge with one of the 
Hubbard bands, a sharp quasi-particle peak emerges at 
w = for sufficiently low T (see inset). 

In a time-of- flight experiment, the two-dimensional 
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— J n-k^dkz/ {2tt), of the three- 
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of the fermions can be measured 0, 17|. Here f{uj) is 
the Fermi function and we have normalized rik such that 
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FIG. 4: (color online) rik (left panel) and njf^ (right panel) 
for k — n/10{nx,ny,nz), nx,y,z = —10, ...,10, plotted as a 
function of ek ~ —2J(cos + cos ky + cosfcz) and e^ = 
— 2 J(cos fci + cosfcy), respectively, for different values of U, 
N and T (upper, middle and lower panels, respectively). For 
a given value of tk a range of values of n(ek) exists. Neverthe- 
less, each curve collapses to a single line in good approxima- 
Inset: An = no,o,o — n^,TT,7T (and An*°^ 
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as a function of U. An and An*°' are largest for the predom- 
inantly metallic phases and smallest for phases with a large 
band- (small U) or Mott-insulating (large U) regions. 
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Fig.[3]shows . For a quantitative analysis of the re- 
sults and for comparison to (future) experiments, we sug- 
gest to plot n^°^ as a function of the tight-binding disper- 
sion e^'' = — 2J(cosfca; -I- cosfcj,), see Fig.^l Surprisingly, 
each curve collapses to good approximation to a single 
line, despite the fact that for a given a range of n^°^ 
exists (as is noticeable in the small scatter of the curves 
for rik). Here it is useful to remember that within DMFT, 
rtk is only a function of et for the homogeneous system. 
Therefore, if a local density approximation (LDA) were 
exactly valid, the collapse to a single curve would be per- 
fect. Therefore our results suggest, that LDA is a very 
good approximation for the analysis of time-of-flight pic- 
tures (but not for other quantities, see below). 

It is an interesting but difficult question whether this 
effect is partially an artifact of DMFT which neglects the 
momentum dependence of the self-energy. While for the 
experimentally relevant temperature range this is proba- 
bly a very good approximation, it is expected to fail very 
close to the metal-insulator transitions at low T. 

Qualitatively, the results of Fig. |4] reflect that local- 
ized electrons in the band- or Mott insulating phase are 
characterized by a momentum-independent rtk while in 
the homogeneous metallic phase displays a jump at 
the Fermi momentum with a height given by the quasi- 
particle weight Z. As the effective local Fermi mo- 
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FIG. 5: (color online) Local spectral function of two non- 
equivalent lattice points at distance 7.794 from the origin at 
the boundary of the Mott insulating region (T — 0.0014Z), 
U = 4.5D). Here the LDA (dot-dashed line) fails completely 
to describe the spectral function but still reproduces the occu- 
pation with high accuracy (?ii/2(9,9,9) = 0.481, ni/2{ii,ii,i) ~ 
0.482, riLDA = 0.496). 



valid close to a sharp domain boundary. Indeed, Fig. [5] 
shows that the LDA fails completely to describe the low- 
energy excitation spectrum at the boundary of the Mott 
insulating region. The coherence peak at the Fermi en- 
ergy arises due to the penetration of the metallic phase 
into the Mott insulator via the Kondo effect. 

For the future it will be interesting to investigate the 
effects of magnetism. In systems with a population im- 
balance we expect that the majority spin will accumulate 
in the Mott insulating regions. These effects will be stud- 
ied in a forthcoming publication. 
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T. Micklitz, D. van Oosten, H. Tjeng, M. Vojta, super- 
computer support by the John von Neumann institute for 
Computing (Jiilich) and the Regional Computing Center 
Cologne and financial support by the SFB 608 of the 
DFG. 



mentum varies smoothly within the trap, all jumps are 
smeared out. However, the slope of n(ek) or the differ- 
ence An*"^ = riQ°Q — 71^°"^ is still a good measure of how 
metallic or insulating the system is. The inset of Fig. 0] 
describes the evolution from a mainly band-insulating 
via a dominately metallic to a Mott insulating regime 
when U/D is increased. Similarly, the middle panel of 
Fig. 3] shows how Ari*°^ increases when at large U/D the 
number of particles and therefore the size of the Mott 
insulating region is reduced (compare with upper inset 
of Fig. [T]). For increasing temperature (lower panel of 
Fig. 01) the destruction of quantum coherence leads to 
a flattening of . Note that is more sensitive to 
changes of T compared to (n^), see Fig. [T] 

Conclusions: In this paper we investigated the sig- 
natures of the Mott transition of fermions in an opti- 
cal trap using a local approximation to the self energy 
(space-resolved DMFT-I-NRG) which allows to treat sev- 
eral thousand atoms. The clearest signature of a Mott 
phase is a plateau in the density profile n{r) of the atoms, 
see Fig. [1] These plateaus are, however, washed out if 
only the column density, J dzn{r), which can be mea- 
sured directly, is considered (not shown). The Mott 
transition is more difficult to observe in a time-of-flight 
experiment. However, the insets of Fig. [?] show that a 
characteristic flattening of n^°^ can be seen when a large 
fraction of the trap becomes a band or Mott insulator. 

Our calculations did not rely on a local density approx- 
imation (LDA) which allows us to investigate whether 
this widely used approximation is valid in the present 
context. It turns out that both density profiles and TOF 
experiments are rather well described by LDA. As dis- 
cussed in the introduction, LDA is expected not to be 
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